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Abstract 



We study four different approximations for finding the profile of discrete solitons in the 
one-dimensional Discrete Nonlinear Schrodinger (DNLS) Equation. Three of them are 
^ • discrete approximations (namely, a variational approach, an approximation to horne- 

ts , clinic orbits and a Green- function approach), and the other one is a quasi-continuum 

approximation. All the results are compared with numerical computations. 

OO 

CN ■ 1 Introduction 
r — 

Since the 1960's, a large number of works has focused on the properties of solitons in the 
Nonlinear Schrodinger (NLS) Equation pp. As it is well known, the one-dimensional NLS 
equation is integrable. Two of the most important discretizations of this equation admit 
discrete solitons. One of these discretizations is known as the Ablowitz-Ladik equation 
[2], which is also integrable. On the contrary, the other important discretization, known 
as the Discrete Nonlinear Schrodinger (DNLS) equation, is not integrable, and discrete 
soliton solutions must be calculated numerically. The DNLS equation has many interesting 
mathematical properties and physical applications [3]. The DNLS equation models, among 
others, an array of nonlinear-optical waveguides [4], that was originally implemented in 
an experiment as a set of parallel ribs made of a semiconductor material (AlGaAs) and 
mounted on a common substrate [5]. It was predicted [6] that the DNLS equation may also 
serve as a model for Bose-Einstein condensates (BECs) trapped in a strong optical lattice, 
which was confirmed by experiments [7i]. In addition to the direct physical realizations in 
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terms of nonlinear optics and BECs, the DNLS equation appears as an envelope equation 
for a large class of nonlinear lattices (for references, see [9], Section 2.4). Accordingly, the 
solitons known in the DNLS equation represent intrinsic localized modes investigated in 
such chains experimentally [10] and theoretically [11[ I12j. In this context, previous formal 
derivations of the DNLS equation have been mathematically justified for small amplitude 
time-periodic solutions in references |13j . 

In this paper we will consider fundamental solitons, which are of two types: Sievers- 
Takeno (ST) modes, which are site-centered [2], and Page (P) modes, which are bond- 
centered [15] (see also Fig. [[]). They can also be seen, respectively, as discrete solitons 
with a single excited site, or two adjacent excited site with the same amplitude. The 
DNLS equation is given by 

iu n + e (u n+ i + u n _i - 2u n ) + j\u n \ 2 u n = 0, (1.1) 

where u n (t) are the lattice dynamical variables, the overdot stands for the time derivative, 
e > is the lattice coupling constant and 7 a nonlinear parameter. We look for solutions 
of frequency A having the form u n (t) = e %ts± v n . Their envelope v n satisfies 



Throughout this paper, we assume 75 > and choose 7 = e = 1 without loss of generality, 
as Eq. (II. 2D can be rescaled. We also look for unstaggered solutions, for which, A > 
(staggered solutions with A < can be mapped to the former upon a suitable staggering 
transformation v n = (—l) n v n ). Furthermore, we restrict to real solutions of (|1.2j) . which 
yield (up to multiplication by expiO) all the homoclinic solutions of (jl.2p [16J. Homoclinic 
solutions of (II. 2h can be found numerically using methods based on the anti-continuous 
limit [11] and have been studied in detail (first of all, in one-dimensional models, but many 
results have been also obtained for two- and three-dimensional DNLS lattices) [3]. 

The aim of this paper is to compare four different analytical approximations of the 
profiles of ST- and P-modes together with the exact numerical solutions. These analytical 
approximations are of four types: one of variational kind, another one based on a poly- 
nomial approximation of stable and unstable manifolds for the DNLS map, another one 
based on a Green-function method, and, finally, a quasi-continuum approach. 

2 Discrete approximations 
2.1 The variational approximation 

Equation (11.21) can be derived as the Euler-Lagrange equation for the Lagrangian 



The VA for fundamental discrete solutions, elaborated in Ref. [T7] (see also Ref. [18]) was 
based on the simple exponential ansatz , 



Av n + e (v n+1 + v n -i - 2v n ) + ^\v n \ 2 v n = 0. 



(1.2) 
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Figure 1. Discrete soliton profiles with A = e = 7 = 1. Left panel corresponds to a ST-mode, 
and right panel, to a P-mode. 




Figure 2. Dependence, for ST-modes, of v (left panel) and v\ (right panel) with respect to A. Full 
lines correspond to the exact numerical solution and dashed lines to the variational approximation. 




Figure 3. Dependence, for P-modes, of v n (left panel) and v\ (right panel) with respect to A. Full 
lines correspond to the exact numerical solution while dashed lines correspond to the variational 
approximation. 
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where v^ T denotes ST-modes, while is for P-modes, with variational parameters A\, 
A2, a\ and a 2 (which determine the amplitude and inverse size of the soliton). Then, 
substituting the ansatz in the Lagrangian, one can perform the summation explicitly, 
which yields the effective Lagrangian, 

rST »r/r.i a on A? tanh 2 ai tP k , ( 2(1- cosh a 2 ) . \ A" 2 2 , 
= M(2sech fll -A-2)+ 2 \ aQh2ai , = A" 2 ^-L___Z_ - Aj +^ tanha 2 

(2.3) 

The norm of the ansatz (|2.2p . which appears in Eq. (|2.3p . is given by Af = J2n=-oo v n- ^ n 
particular, for the ST- and P-modes, 

Mi = A^cothai, A2 = A|/sinha 2 . (2.4) 

The Lagrangian (12, 3p gives rise to the variational equations, OL^T /dMi = dL^S / da\ = 
0, and dL^ s /dAf2 = dL^/da^ = 0, which constitute the basis of the VA [19]. These 
predict relations between the norm, frequency, and width of the discrete solitons within 
the framework of the VA, namely 

4 cosh a\ sinh 2 2a\ . 8(1 — cosh a 2 + sinha 2 ) cosh 2 a 2 

Ai = . , , r-r-^ — , A 2 = — : (2.5) 

sinn 4a 1 — smn 2ai sinn a 2 + cosn a 2 

. r tanh 2 ai , 2(1 — cosh a 2 ) 1,, , , 

A = 2 sech ai - 1 )+A V , n , A = -r\ ^- + -M 2 tanh a 2 . 2.6 

tanh 2ai smn a 2 + cosn a 2 2 

These analytical predictions, implicitly relating J\f and A through their parametric de- 
pendence on the inverse width parameter a, will be compared with numerical findings 
below. In Figs. [2] and [3l we compare the approximate and exact values of the highest 
amplitude site and the second-highest amplitude sites (i.e. vq and v\, which can be eas- 
ily calculated from (|2.5p once N and a are known) with respect to A for both ST- and 
P-modes. We can observe that the variational approach captures the exact asymptotic be- 
havior as A — > +00. Indeed as a\ — > +00 in approximation (|2.2p one obtains A ~ Ai ~ e ai 
and A\ ~ VAi ~ \/A. Thus Vq T ~ y/A as A — > +00 which is indeed the asymptotic be- 
havior of the exact ST-mode. On the contrary, the variational approximation errs by a 
small multiplicative factor (-^= ~ 1.1) as A — ► (i.e., effectively approaching the contin- 
uum limit). This can be seen taking the limit a\ — > in approximation (|2.2p . One has 
M\ ~ 8ai, A ~ —a 2 + f Ai ~ 3a 2 and A\ ~ 2\f2a\ ~ ^=\/2A, while the amplitude of the 

continuum hyperbolic secant soliton of the integrable NLS is A = \/2A [see also below]. 
Notice that the P-mode also has the same A — * limit (and therefore errs by the same 
factor) . 

2.2 The homoclinic orbit approximation 
2.2.1 The DNLS map 

The difference equation (jl.2p can be recast as a two-dimensional real map by defining 
y n = v n and x n = v n . l (23 EH M, M, M\ : 



Vn+i = -Vn + (A + 2)y n - x, 



3 , ,A , ON.. _ (2-7) 
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Figure 4. Homoclinic tangles for A = 0.4, A = 0.6. 



For A > 0, the origin x n = y n = is hyperbolic and a saddle point, which is checked upon 
linearization of the map around this point. Consequently, there exists a 1-d stable and 
a 1-d unstable manifolds emanating from the origin in two directions given by y = X±x, 
with 



(2 + A) ± ^A(A + 4) 



(2. 



The eigenvalues A± satisfy A 2 — (A + 2) A + 1 = and A + = Al 1 > 1. The stable and 
unstable manifolds are invariant under inversion as it is the case for eq. (|2.7|) . Moreover, 
they are exchanged by the symmetry (x, y) >— » (y, x) (this is due to the fact that the map 
(|2.7|) is reversible; see e.g. [IB] for more details). Due to the non-integrability of the DNLS 
equation, these manifolds intersect in general transversally, yielding the existence of an 
infinity of homoclinic orbits (see Figs. 0]and[5]). Each of their intersections corresponds to a 
localized solution, which can be a fundamental soliton or a multi-peaked one. Fundamental 
solitons, the solutions we are interested in, correspond to the primary intersections points, 
i.e. those emanating from the first homoclinic windings. Each intersection point defines 
an initial condition (xo,yo), that is, (v-i,Vq), and the rest of the points composing the 
soliton are determined by application of the map. 



2.2.2 The polynomial approximation to the unstable manifold 

The first windings of the stable and unstable manifolds can be approximated by third order 
polynomials. Actually, only one of them is necessary to be determined, as the other one 
is determined taking into account the symmetry x «-> y. We proceed then to approximate 
the local unstable manifold Wj" (0). Taking into account its invariance under inversion, 
it can be locally written as a graph y = f(x) = Xx — ax 3 + 0(\x\ 5 ) with A = A + given by 
(|2.8p . For x ~ 0, the image of (x, f(x)) under the map (|2.7p also belongs to Wj" c (0), thus 
-f(xf + (A+2)f(x)-x = f(f(x))Vx « 0. This yields [A 3 +a(A+2-A-A 3 )]x 3 +0(|x| 5 ) = 
0, Vx ss 0. Hence a = -A 3 /(A + 2 - A - A 3 ) = A 4 /(A 4 - 1). The local unstable manifold 
is approximated at order 3 by 

W u : y = Xx - ^—^x\ (2.9) 
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A=1 A=3 




Figure 5. Homoclinic tangles for A = 1 and A = 3. 




Figure 6. Numerical exact unstable manifold (full line) and its approximation by Eq. (|2.9| 
(dashed line) for A = 1 (left panel) and A = 3 (right panel). The fit is so accurate in the latter 
that both curves are superimposed. 

and, by symmetry, the stable manifold is approximated by: 

W s : x = Xy - j^jV 3 - (2.10) 

In Fig. [U the numerical and approximated unstable manifolds for A = 1 and A = 3 are 
compared. It can be observed that the fit is better when A increases. The approximation 
breaks down for small A because the origin is not a hyperbolic fixed point for A = 0. 

2.2.3 Approximate solutions via approximate invariant manifolds 

Once an analytical form of the unstable and stable manifold is found, discrete solitons 
profiles (or, concretely, vq and i>-i) can be determined as the intersection of both manifolds. 
The polynomial form of (j2.9[) is not sufficient in practice to obtain good approximations of 
the whole soliton profile, due to sensitivity under initial conditions. However, it provides 
a good approximation near the soliton center. Some intersections of W s and W u can be 
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Figure 7. (Left panel) Approximated stable and unstable manifolds for A — 2 showing the main 
intersections. (Right panel) Pitchfork bifurcation arising in the homoclinic approximation when A 
is varied. ST-modes (full lines) bifurcate with the P-mode (dashed line) at A = 0.5. 




A A 



Figure 8. Same as Fig. [2] but with dashed lines corresponding to approximation (|2.12p . 



approximated by: 

x = a(a.-^)-^(a,-_^)\ (2 . n) 

This equation has nine solutions (see Fig. [7^). One of them {x = 0), corresponds to the 
origin. Once this solution is eliminated, the reminder equation is a bi-quartic one. Thus, 
if x = £ is a solution of (|2.1ip . x = — £ is also a solution: this is due to the fact that ±v n 
is a solution of (|1.2p . Solutions x = £i, x = £21 % = £o and x = £3 i n Fig. [7] correspond to 
the positive solutions of (|2.1ip . The point x = £0 is in the bisectrix of the first quadrant 
and corresponds to the P-mode (i.e. Vq = £o)> an d the point x = £3 lies in the bisectrix of 
the fourth quadrant and corresponds to a twisted mode (i.e. a discrete soliton with two 
adjacent excited sites with the same amplitude and opposite sign). Setting y(£o) = £0 and 
2/(6) = -6 m (ZD), one obtains £0 = A~V(A - 1)(A 4 - 1), £3 = A" 2 y/{\ + l)(X 4 - 1). 
Upon elimination of the roots x = £0 and x = £3 from (|2.1ip . £1 and £2 can be calculated 



8 



J Cuevas, G James, P G Kevrekidis, B A Malomed and B Sanchez-Rey 




Figure 9. Same as Fig. [3] but with dashed lines corresponding to approximation (|2. 13|) . 
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Figure 10. Same as Fig. [2] but with dashed lines corresponding to approximation (|2.18|) . 

as solutions of a quadratic equation. Thus, 

Ci = A-V(A4-l)(A- V / A2-4)/2, 6 = A- 2 V / (A4-l)(A + v /A 2 -4)/2. (2.12) 

These solutions are related with the ST-mode as Vq T = £2 and vf T = £1. On the other 
hand, for the P-mode, Vq = Co, and, vf should be determined by application of the map 
(12771) . This yields 

^o P = A-V(A-l)(A 4 -l), uf = A~ 6 (A 3 + A - l)v/(A - 1)(A 4 - 1). (2.13) 

In Figs. [8]and[9l the values of vq and v\ obtained through the homoclinic approximation 
are represented versus A and compared with the exact numerical results. It can be observed 
that, for ST-modes, no approximate solutions exist for A < 0.5. For A = 1/2 (i.e. A = 2), 
the points (£1,^2) and (£2,^1) disappear via a pitchfork bifurcation at (£o>£o) ( see Fig- 
[7b). This artifact is a by-product of the decreasing accuracy of our approximations as 
A — > 0; as discussed before, the ST-mode should exist for all values of A > 0. 

2.3 The Sievers-Takeno approximation 

A method to approximate solutions of (jl.2p has been introduced by Sievers and Takeno, 
for a recurrence relation similar to it but with slightly different nonlinear terms |14j . 
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This approach has been generalized to the d-dimensional DNLS equation in reference 
|23j. In what follows we briefly describe the method, incorporating some precisions and 
simplifications. Setting v n = Vof] n , equation (|1.2p becomes 

%+i - 2ri n + rj n -i = Kr} n - vqT}\, (2.14) 

with = ry„, r/g = 1. Setting n = in (|2.14p we obtain in particular 

v 2 = A + 2(l- m ). (2.15) 

Equation (|2.14p can be rewritten as a suitable nonlocal equation using a lattice Green 
function in conjunction with the reflectional symmetry of r\ n and equation f|2.15j> . This 
yields for all n > 1 

r) n = [A + 2(1 - m ) + E ^ ( A ~ |n ~ fe| + A ~ n ~*)> ( 2 - 16 ) 

fe>i 

where A = A+ is given by (]2.8p . Problem (|2.16p can be seen as a fixed point equation 
{r/} = F\({r]}) in (N*). Noting B e the ball ||{^}||^ 00 (n*) — e J the map Fa is a contraction 
on B e provided e is sufficiently small and A is greater than some constant Ao(e). In that 
case, the solution of (|2.16p is unique in B e by virtue of the contraction mapping theorem 
and it can be computed iteratively. Choosing {77} = as an initial condition, we obtain 
the approximate solution 

Vn^(F A (0)) n = ^L\- n , n>l. (2.17) 

Obviously the quality of the approximation would increase with further iterations of -Fa- 
Using (|2.17p and (|2.15p in the limit when A is large, we obtain 

v n « (A + 2) 1 / 2 A- |n| (2.18) 

since A ~ A as A ^ +00. The values of vq and v\ in this approximation are compared 
with the exact numerical results in Fig. [TUJ We observe that the approximation captures 
the asymptotic behaviour of vq and v\ for A — ► 00. 



3 The quasi- continuum approximation 

As it can be concluded from previous sections, none of the established approximations 
perform well for A close to zero (although the VA is notably more accurate than the 
invariant manifold and Sievers-Takeno approximation) . A quasi-continuum approximation 
could be used to fill this gap. To this end, we follow Eqs. (13) and (14) of Ref. [M]- Then 
the ST- and P-modes can be approximated by the continuum soliton based expressions: 

= \/2Asech (nVX), vl = \/2Asech [(\n + 1/2| - 1/2)VA]. (3.1) 

These expressions lead to the results shown in Figs. [TT]and[T2j Naturally, this approach 
captures the asymptotic limit vq ~ \/2A when A — > 0, but fails increasingly as A grows. 
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Figure 11. Same as Fig. [5] but with dashed lines corresponding to approximation (|3.ip . 




Figure 12. Same as Fig. [3] but with dashed lines corresponding to approximation (|3.1D . 

4 Summary and conclusions 

In Figs. [13] and O the results of the paper are summarized. To this end, a variable, giving 
the relative error at site n, is defined as: 

Rn = log 10 |« PPr ° X " C aCt )/C aCt | • (4-1) 

We can generally conclude that the variational approximation offers the most accurate 
representation of the amplitude amplitude of the Page mode v£ at the two sites n = 
and n = 1 with some small exceptions. These involve some particular intervals of A where 
the homoclinic approximation may be better and also the interval sufficiently close to 
the continuum limit, where the best approximation is given by the discretization of the 
continuum solution. Similar features are observed for the approximation of the Sievers- 
Takeno mode at site n = 0. However, a different scenario occurs for this mode at site 
n = 1, since the homoclinic approximation gives the best result for A > 1.5. As A goes to 0, 
the Sievers-Takeno, variational and quasi-continuum approximations give successively the 
best results in small windows of the parameter A. Notice that in the interval A € (0, 0.5] 
neither the variational, nor the homoclinic approximation are entirely satisfactory. The 
latter suffers, among other things, the serious problem of producing a spurious bifurcation 
of two ST modes with a P-mode. On the other hand, for larger values of A (i.e., for 
A > 0.5), the quasi-continuum approach is the one that fails increasingly becoming rather 
unsatisfactory, while the discrete approaches are considerably more accurate, especially for 
A > 2, when their relative error drops below 1% (with the exception of the Sievers-Takeno 
approximation of Vq T , which only reaches this precision for A > 10). 
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A A 

Figure 13. Representation of variable R denned in (|4.1j) versus A for ST-modes. Full lines corre- 
spond to the variational approach; the dashed line corresponds to the homoclinic approximation; 
the dash-dotted lines to the continuum approximation; and the dotted line to the Sievers-Takeno 
approximation. 




A A 

Figure 14. Representation of variable R defined in (|4.1| versus A for P-modes. Full lines cor- 
responds to variational approach; dashed line, to the homoclinic approximation; and dash-dotted 
lines, to the continuum approximation. 



We hope that these results can be used as a guide for developing sufficiently accurate 
analytical predictions in different parametric regimes for such systems. It would naturally 
be of interest to extend the present considerations to higher dimensions. However, it 
should be acknowledged that in the latter setting the variational approach would extend 
rather straightforwardly, while the homoclinic approximation is restricted to one space 
dimension and the other approximations would become more technical. 
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